AD - A 1 4 1  731  ANALYSIS  OF  AN  INTEGRATION  METHOD  FOR  RAPIDLY  1  /  1 

OSCILLATING  INTEGRANDS! U )  WISCONSIN  UNIV-MADISON 
MATHEMATICS  RESEARCH  CENTER  D  LEVIN  ET  AL .  APR  84 
UNCLASS  I F I  ED  MRC-TSR-?fi70  DAAG7q-Rn-C  0041 _ F/R  19/1  Ml 


Mathematics  Research  Center 
University  of  Wisconsin— Madison 
610  Walnut  Street 
Madison,  Wisconsin  53705 


'^April  1984 

Q_ 

o 

CJ> 


LxJ 


(Received  July  18, 


1983) 


Sponsored  by 

U.  S.  Army  Research  Office 
P.  O.  Box  12211 
Research  Triangle  Park 
North  Carolina  27709 


Approved  for  public  release 
Distribution  unlimited  L'  * 


ELECTEp 

JUU  0  1  1984 


84  05 


31 


064 


E 


UNIVERSITY  OF  WISCONSIN-MADISON 
MATHEMATICS  RESEARCH  CENTER 


ANALYSIS  OF  AN  INTEGRATION  METHOD  FOR  RAPIDLY  OSCILLATING  INTEGRANDS 
D.  Levin*,  L.  Reichel  and  C.  Ringhofer 


Technical  Summary  Report  #2670 
April  1984 


v  ABSTRACT 

"SA  collocation  method  for  evaluating  integrals  of  rapidly  oscillating 
functions  is  analysed.  The  method  is  based  on  the  approximation  of  the 
antiderivative  of  the  integrand  as  the  solution  of  an  ordinary  differential 
equation  via  polynomial  collocation.  Numerical  examples  are  presented. 
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SIGNIFICANCE  AND  EXPLANATION 


Special  integration  methods  have  to  be  used  for  the  integration  of  highly 
oscillatory  functions.  In  the  present  paper  we  describe  such  a  method  and 
determine  the  order  of  convergence. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
suasiary  lies  with  MFC,  and  not  with  the  authors  of  this  report. 


ANALYSIS  OF  AN  INTEGRATION  METHOD  FOR  RAPIDLY  OSCILLATING  INTEGRANDS 


D.  Levin*,  L.  Reichel  and  C.  Ringhofer 

1.  INTRODUCTION 

In  (3]  a  collocation  method  for  computing  integrals  with  rapidly  oscillatory 
integrands  was  presented,  and  numerical  examples  Indicated  fast  convergence.  In  the 
present  paper,  we  analyze  the  method,  show  convergence  rates,  and  discuss  the  choice  of 
collocation  points. 

we  are  to  compute  integrals  at  the  form 

(1.1)  I  -  ]  f(x)elq<X,/Edx  , 

a 

where  -«•  <  a  <  b  <  •,  i  •  /-T,  q(x)  is  a  real-valued  function,  f(x)  is  a  complex¬ 
valued  function  and  e  is  a  constant  0  <  e  *  b  -  a.  He  assume  that  f(x)  and  q(x) 
vary  moderately  on  [a,bl .  Specifically,  we  assume  there  are  constants  5 1 ,S  2  such  that 
0  <  $i  <  <x)|  <  and  >e  holds.  The  integration  method  is  based  on  the 

observation  that  if  a  function  p(x)  satisfies 

(1.2)  (p(x)«iq(x,/E)  -  f(x)elqtx)/E,  a  <  x  t  b  , 
then  the  integral  is  given  by 

(1.3)  I  -  p(b)elq,b,/E  -  p(a)eiqU,/E  . 

From  (1.2)  it  follows  that  p(x)  is  any  solution  of 

(1.4)  p'  <x)  +  q'  (x)p(x)  -  f(x) ,  a  <  x  <  b  . 

Equation  (1.4),  we  solve  by  approximating  p(x)  by  piecewise  polynomials,  which  we 
determine  by  collocation.  Details  ae  well  as  convergence  results  are  provided  in  Section 

2.  Section  3  contains  numerical  examples. 


•School  of  Mathematical  Sciences,  Tel-Aviv  University,  Raraat-Aviv,  Israel 


Sponsored  by  the  united  States  Army  under  Contract  No.  DAAG2R-80-C-0041 . 


4 


2.  CONVERGENCE  PROPERTIES 

Lot  a  ■  z1  <  z2  <  •••  <  *m+i  “  b  h*  a  partition  of  Ca,b] ,  and  for  notational 
convenience  assume  that  all  z j  are  equidistant  with 

(2.1)  h  -  *j+1  “  *j.  j  “  • 

and  h61  >  e.  Cn  each  subinterval  z^) ,  we  approximate  p(x)  by  a  k-th  degree 

polynomial,  which  we  determine  by  collocation  in  k  +  1  points.  This  gives  ua  a 

2  k-1 

convergence  order  of  0(£  h  ),  and  as  example  2.1  below  shows,  this  result  is  sharp. 
From  lemma  2.1  below  it  follows  that  jl|  -  0(e),  why  the  relative  accuracy  of  the  method 
is  0<ehk  ) .  We  next  describe  the  collocation  method.  Choose  k  +  1  reference  points 
w ^ ,  0  «  w0  <  w,  <  ...  <  -  1.  The  collocation  points  w^  in  [Zji.1,Zj],  we  obtain  by 

the  linear  map 

(2.2)  ws  ”  zj  +  hws'  8  “  0(1)k,  j  “  1(1)m  . 

Let  i  ,(x)  be  the  Lagrange  polynomials  defined  by  w-L  s  “  0(1)k,  i.e. 

"3  ■ 


(2.3) 


ln,j(x> 


k  X  -  MJ 

_  _  8 

n  .  .  , 

s-0  wJ  —  wJ 

.  n  s 

s*n 


n  -  0(i)k  . 


Determine  the  polynomial 

(2.4) 


PjU) 


k 

l 

n-0 


to  “"3l»3(X) 


from 

(2.5)  p’(x)  +  ^  q'(x)p^(x)  -  f(x),  X  -  wj*,w^,...,wj|  . 

The  computed  p^(x)  approximates  p(x)  on  [*j,Zj+1].  Analogously  to  (1.3) 

iq(z  )/e  iq(z  )/e 

(2.6)  Ij  s*  Pj{*j+l,e  ^  -  Pj(*j,e 


is  an  approximation  of 

(2.7) 

approximates  I. 


j  j+1  f(x)elq(x)/edx.  and 

*3 

jk  1'  xk 
i-1  i 


-2- 


Theorem  2. 1 


Let  I,  lK,  h,  Pj  etc.  be  defined  by  (1.1)  and  (2.1)-(2.7).  Then 

!  I  -  I  1  <  ce  h  , 

where  c  is  a  constant  independent  at  e  and  h. 

To  proove  theorem  (2.1)  we  first  need  some  additional  lemmast  As  we  will  see  it  is 
essential  for  our  method  that  there  exists  a  sufficiently  smooth  solution  of  (2.5),  i.e. 
that  there  exists  a  solution  which  has  sufficiently  many  derivatives  bounded  uniformly  as 
E  goes  to  zero.  This  is  shown  in 


Lemma  2. 1 

2*-l  2k 

Let  f(x)  e  C  [Zj,z  and  q(x)  e  C  [z^.z^].  Then  there  is  a  solution  v(x) 

(2.8)  v'(x)  +  j  q' (x)v(x)  »  f(x),  z^  <  x  <  z^+1  , 
which  satisfies 

.a 

(2.9)  | — ^  (x)|  <  ce  on  »  “  0(1)x 

dx 

where  c  is  a  constant  independent  of  e . 

Proof.  Let 

x 

vk(x)  l  v  (x)E8  , 
a«1 


v^x) 


f(x) 

'  iq'(x)' 


v  (x) 

8 


v;-i(*> 

lq’(x)' 


2(1)X  . 


Then  v^  satisfies 

v;(x)  +  ■-  q'(x)vK<x>  -  f(x,  +  EkvJ(x) 

and 


Now 


dx 


B  (X)  -  0(E),  Z^  <  X  < 


S  -  0( 1 )K  . 


v ( x )  VK(Zj 


).iq,aj)/'.-iq(X>A  +  e-iq(x)/e  ]*  eiq<t)/ef(t, 


dt 


-3- 


solves  (2.8).  Introduce  e(x)  :•  v(x)  -  V^lx).  Then 


e*(x)  +  ~  q*  (x)e(x)  -  v^(x)ek,  ed^)  -  0  . 


Then  e(x)  can  be  written  as 


e(x)  -  eViq,x,/e  I  ei<j(t,Av'(t)dt  . 


J  eiq<t,/Ev^(t)dt  -  0(e),  Sj  <  x  <  r^+1  , 


It  follows  that  e(x)  -  0(e  )  on  Repeated  differentiation  of  (2.11)  shows 

that  d  8  (x)  is  at  most  0(e'C+1  fl),  s  »  0(1)x.  By  construction  V  (x)e  1  has  < 

.  a  x 

dx 

derivatives  bounded  independently  of  e.  The  lemma  follows  from  v(x)  *  V  (x)  +  e(x). 
Remark  2. 1 

The  smoothness  requirements  on  f(x)  and  q(x)  are  only  sufficient.  It  is  simple  to 
construct  functions  where  v(x)  has  more  continuous  derivatives  than  f(x)  and  g(x). 

See  section  3. 

Next  we  estimate  the  error  in  Ik  (see  (2.6))  on  each  subinterval  tzj'zj  +  1^*  I*t 
v(x)  solve  (2.8).  Then 


f(x)eiq(x,/Edx  -  Ik| 


l,v<8j+1*  “  p)(Vi,)# 


iq(*J+i)/e 


iq(  Z  )/t 

(vd  )  -  d^ )  )•  3  |  < 


<  |vdj+,)  -Pj<zj+1)|  +  |v<zj)  -  Pj  <*j )  I  • 

Hence,  it  suffices  to  bound  v(x)  -  Pj(x)  of  z ^  and  zj+i* 
liSfa  2.  2 

Let  v(x)  and  Pj(x)  have  the  same  meaning  as  in  (2.11).  At  the  collocation 
points  wi,  s  ■  0(1)k,  in  fzj'zj+lJ'  (see  (2.2)) 


i  1 1  llliT  il  I  III  I  i 


-4- 


(2.12) 


|v(w^)  -  p.(w^)|  <  ec  max  | r  |  , 
*  3  ■  <Xs<k  * 


holds  whsra  c  Is  a  constant  Independent  of  e  and  h,  and  the  r#  are  defined  as 

*4  4 

follows!  let  Pj(x)  be  the  polynomial  of  degree  <  k  satisfying  p^(w')  -  v(w'), 
s  ■  0(1  )k.  We  define  r(  as  the  residuals 

(2.13)  i-  v'(w^)  -  p*(w^),  ■  “  0(1)k  . 

Proof 

With  the  representation  (2.4)  of  p^tx),  the  collocation  equations  (2.S)  become 

(2.14)  (L  +  £  E»a  -  t  , 

where 

D  I  -  diagdq'  (w^)  ,iq'  (w^) , .  ..,iq*  (w^) ) 


'Vo*  *  *  ‘ 


L  i- 


t'  (wj)  ♦  •  •  t’Jw;’) 


Oj  k 


kj  k 


f  (f(w^ ),..., f(w£))T 


a  »■  (i 


OJ' 


'V 


4  IT 

Let  v(x)  be  any  solution  of  (2.8),  and  let  v  i-  (v(w' ) , . . . ,v(w') )  .  We  write 
P'(wJ)  +  *1’  <w*)Pj  <w£ )  -  «wj)  +  rg,  s  •  0 ( 1  )k,  as 


(2.15) 


(L  +  —  n)v  -  _f  +  £  , 


T  lii 

where  r  «  (r(J/r1,...,rk)  .  Since  the  elements  of  L  are  O(-),  |q'  |  >  and  e  <  6,h 

the  matrix  L  +  ^  d  is  invertible,  and 


v  -  a  »  c(el  4-  D)  'r  -  ed_1  \  (-LD  ’e)8^  . 

s»0 


-1 


This  proves  the  1« 


-5- 


The  next  lemma  bounds  r.  The  result  Is  Mil-known  end  a  simple  proof  by  Itolle's 
theorem  is  omitted . 

Lemma  2.3 

Let  p(x)  be  the  polynomial  of  degree  <  k  that  interpolates  g(x)  C  C*t+1[rj,Sj+1) 
at  points  *  •  •  •  Then  for  any  x  e  f*j<*j+il» 


|p'(x)  -  g*  (x)|  <  chk  max 


.k+1  ,  . 

id  I 


dx 


k+1 


where  c  is  a  constant  independent  of  h. 

We  are  now  in  a  position  to  give  an  estimate  for  the  integration  error  for  one 
subinterval . 

Theorem  2.2 

Let  Ik  be  defined  by  (2.6),  and  assume  that  f(x)  and  g(x)  satisfy  the  smoothness 
requirements  of  lemma  2.1.  Then 


z 

(2.16)  |l*  -  ]  f(x)eiq<x>/edx|  <  ce2hk  , 

where  o  is  a  constant  Independent  of  e  and  h. 

Proof.  Taking  the  smooth  solution  of  lemma  2.1  with  x  •  k  +  1,  its  k  +  1st  derivative 

in  continuous  and  0(e).  By  lesna  2.3,  rs  introduced  in  lemma  2.2  satisfies 

max  Ir  I  <  cehk  for  some  constant  c  independent  of  e  and  h.  Therefore  the  theorem 
0<s<k  * 

follows  by  (2.11)  and  (2.12). 

The  proof  of  theorem  2.1  is  now  trivial.  It  follows  from  theorem  2.2  and  summation 
over  all  m  -  (b  -  a)/h  subintervals.  The  following  example  shows  that  the  error  estimate 
of  theorem  2.1  is  sharp  if  e  <  hi 1 . 

Kx.  2.1  consider  the  integral 

I  -  J  1  exeix/edx  . 

0 

Divide  [0,1]  into  m  subintervals  [ZjfSj+f]  of  length  h  1/m,  i.e.  ij  “  hj. 


-6“ 


j  -  0(1 )b.  For  future  reference,  we  write  I  as 


(2.17) 


m*1  *j+1  ,  „  -  .  <w,  *-1  z  iz  /e 

l  ]  e*ei*/  dx  -  -j-~  (eheih/e  -  1)  \  -  3-  j 


j-0  z 


j 


j-0 


We  apply  our  integration  method  with  k  -  1  on  each  subinterval  ^cj,zj+1^‘  We  determine 
a  linear  function 

S  -  X  X  -  z 

p  ( x )  *  d  *  +  d  ■  -  * 

j  oj  h  ij  h 


which  satisfies 


lliis  determines 


pj<V  "I  pj<V  “e8'  s  “  *  1 


m-1  iz^_^,/e  izA/e 

“i  j  "  d0j ' 


X1  ,-T  I-  i+,  \-  -e“J'  , 


j-0 


which  after  some  simplification  can  be  written  as 


(2.10) 


.1  re  h  ih/e  .  _ 

I  »[j-(ee  -1)+e 


2  (,^A 


m-1  z.  iz  /e 
-  1>M  e  3e  3 
j-0 


with  (2.17)  this  gives 


(2.19)  I1  -  I  -  < 


-e _  h  ih/e 


e  e  -  1)  +  e 

II 

m-1  z.  iz./e 


2  .^1  (eih/e  .  ,)]  y  e~3e~""3 


i-l  z4  izye 

j-0 


[0(e3)  +  e^Oth) ]  l  e  'e 

j-0 


The  sum  in  (2.19)  is  bounded  by  0(1/h),  and  therefore 

1 1 1  -  l|  <  0<e2)  +  0(e3/h)  . 

From  (2.19)  it  is  obvious  that  theorem  2.1  gives  a  bound  for  the  decay  of  the  error,  while 
the  error  itself  might  not  tend  to  zero  smoothly  as  h  ♦  0  or  e  ♦  0.  This  is  also 
illustrated  by  the  numerical  examples  of  section  3. 


-7- 


We  conclude  this  section  with  a  suggestion  for  collocation  points.  In  the  proofs  we 


have  only  required  that  the  collocation  points  are  distinct  on  each  subinterval 

and  that  the  end  points  of  each  subinterval  are  collocation  points.  While  the 

|( 

allocation  of  collocation  points  does  not  affect  the  rate  of  convergence  I  ♦  I  as 
h  ♦  0  or  as  e  ♦  0,  the  allocation  does  affect  the  constant  c  in  (2.16)  and  in  theorem 


2.1.  Let  D  denote  the  differentiation  operator,  and  write  (2.8)  symbolically  as 

(1  +  IT-  D)v  “  ifr  f 

or  equivalently,  if  e  is  sufficiently  small  and  f  is  sufficiently  smooth. 


This  shows  that  the  polynomial  Pj(x)  we  determine  by  (2.5)  is  close  to  the  kth  degree 

polynomial  p*(x)  which  interpolates  r— T  f  at  wj,  s  ”  0(1)k,  and  it  suggests  the  us 
j  iq'  8 

of  collocation  points  w^,  which  are  good  for  the  interpolation  problem 
P^(w^)  -  J-  f  (w^ )  (q*  (w^ ) )  1 ,  s  -  0(1)k.  If  we  select  the  wj*  as  the  extended  Chebyshev 

j  S  1  8  8  8 

points 


(2.21) 


.  cos< ( 2t  +  1)x/(2k  +2)),  . 

1>  - cos("x/( 2k  ;  2,) - ]>  t“0(1)k' 


then  p*(x)  will  be  close  to  the  best  polynomial  approximant  to  -j  f(x)(g'(x))  \  see 
de  Boor  (1],  ch.  2  or  Brutman  (2].  The  examples  of  section  3  confirm  that  the  extended 
Chebyshev  points  also  are  good  for  solving  the  collocation  problem  (2.5). 


3.  NUMERICAL  EXAMPLES 

In  all  examples  the  Interval  la  [0,1).  The  ra  eubintervala  [ZyZy^l  are  defined 
by  «  h y  h  -  jj,  J  -  0(1)ra.  k  denotea  the  degree  of  the  polynomial  Pj(x).  If 
nothing  elee  stated,  the  collocation  points  ^/s-O  are  the  extended  Chebyshev  points 
(2.21).  All  computations  were  done  on  a  VAX  11/260  in  double  precision  arithmetic,  i.e. 
with  12  significant  digits. 

Ex.  3.1  Let  q(x)  tm  x  +  x2  and 


f(x) 


i 


1  +  2x 

(x  -  0.S)2  +  0.1 


( (x 


2x  -  1 _ 

0 .5 ) 2  +  0.1)2 


e . 


Then  (1e4)  is  solved  by  p(x) 


e((x  -  0.5)2  +  0.1)"1. 


The  table  shows  the  dependence  of 


h 

k 

e 

l*5  (rounded) 

|i  -  ik| 

1 

5 

n 

I 

o 

-3.9093  ♦  10-3  +  2.6642  •  10'31 

7.3  •  10-( 

1/2 

5 

10-3 

-3.9070  •  10-3  +  2.6573  •  10~3i 

1.4  •  10"1 

1/4 

5 

10"3 

• 

5.0  •  10-1 

1/0 

5 

10-3 

m 

4.7  •  10"’ 

1/16 

5 

10"3 

m 

1.7  •  10" 

1 

5 

IQ"6 

-6.9997  •  10"7  - 

1.8735  •  10'6i 

1.0  • 

10"11 

1/2 

5 

io-« 

m 

1.5  • 

io-12 

1/4 

5 

io-« 

m 

8.8  • 

io-14 

1/8 

5 

10-6 

m 

3.7  • 

10*15 

1/16 

5 

IQ'6 

m 

2.3  • 

IO"16 

1/2 

20 

10”3 

-3.9070  •  10-3  + 

1.6573  •  10*3i 

2.9  • 

10-14 

1/2 

20 

10*3 

-6.9997  •  ID'7  - 

1.8735  •  10”6i 

1.4  • 

io-19 

1/2 

20* 

10*3 

-3.9070  ♦  10"3  + 

2.6573  •  10*3i 

1.3  • 

10-11 

•equidistant 

collocation 

point 

We  note  that  high  accuracy  can  be 

obtained  not  only  by 

decreasing  h, 

but  also  by 

increasing 

k.  the  latter  strategy  can  be  competitive 

since  the  strong  diagonal  dominance 

of  aystesi  (2 

.14)  allows 

iterative 

solution  in  0(k2)  operations. 

Ex.  3.2  Let  g(x)  j-  x  +  x2  and  f(x)  i-  (1+2x)(x - )2sign(x - -)i  -  2(x — )e . 

12  1  /2  /2  /2 

Then  p(x)  i“  (x - )  sign(x - )  solves  (1.4). 

✓  2  /2 


h 

k 

e 

I*  ( rounded ) 

|l  -  Ik| 

1/8 

2 

10“3 

4.6847  •  10-4  +  7.9750  •  10-5i 

3.63  •  10-8 

1/16 

2 

10-3 

4.6849  j  10"4  7.9770  •  10_5i 

1.85  •  10"8 

1/4 

8 

10-3 

4.6848  • 

10'4  r 

/  .9788  • 

io-5i 

4.47  •  10 

1/8 

8 

ID"3 

4.6848  • 

o 

1 

* 

7.9783  • 

io-5i 

2.74  •  10' 

1/4 

8* 

10"3 

4.6849  • 

10"4  + 

7.9809  • 

10“5i 

2.9  •  10_l 

•equidistant  collocation  points • 

The  example  shows  the  robustness  of  the  method. 

The  last  example  shows  that  the  conditions  of  letmoa  2.1  are  not  necessary. 

Ex.  3.3  Let  q(x)  |x - and  define  q’(— r)  ”  0-  Let 
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Then  p(x)  is  the  same  as  in  Example  3.1. 
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